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ABSTRACT. We present algebraic multilevel iteration (AMLI) methods for isogeometric discretization of 
scalar second order elliptic problems. The construction of coarse grid operators and hierarchical comple- 
mentary operators are given. Moreover, for a uniform mesh on a unit interval, the explicit representation of 
B-spline basis functions for a fixed mesh size h is given for p = 2, 3, 4 and for C°- and C p_1 -continuity. 
' J— " ' The presented methods show h- and (almost) p-independent convergence rates. Supporting numerical re- 

sults for convergence factor and iterations count for AMLI cycles (V-, linear W-, and nonlinear W-) are 

f*^ ■ provided. Numerical tests are performed, in two-dimensions on square domain and quarter annulus, and in 

r\] ' three-dimensions on quarter thick ring. 
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1. Introduction 

The IsoGeometric Analysis (IGA), proposed by Hughes et al. in [2 ], has received great attention in 
the computational mechanics community. The concept has the capability of leading to large steps forward 

^*«i ' ■ in computational efficiency since effectively, the process of re-meshing is either eliminated or greatly 

suppressed. The geometry description of the underlying domain is adopted from a Computer Aided De- 
sign (CAD) parametrization which is usually based on Non-Uniform Rational B-splines (NURBS), and 

*£ the same basis functions are employed to approximate the physical solution. Since its introduction, IGA 

techniques have been studied and applied in diverse fields, see e.g., [1, 8, 9, 16, 17, 20, 26, 28]. Moreover, 
some theoretical aspects such as approximation properties, condition number estimates have been stud- 
ied, see [7, 10, 14, 25]. The isogeometric methods, depending on various choices of basis functions, have 
shown several advantages over standard Finite Element Methods (FEM). For example, some common 
geometries arising in engineering and applied sciences, such as circles or ellipses, are exactly repre- 
sented, and complicated geometries are represented more accurately than traditional polynomial based 
approaches. When we compare NURBS based isogeometric analysis with standard Lagrange polynomi- 

O ■ als based finite element analysis, it leads to qualitatively more accurate results [19]. Another limitation 

t^j- \ of finite element analysis is that it suits well for C° continuous interpolation, but for C 1 or higher order 

interpolation finite elements are complicated and expensive to construct. IGA offers C p_fc -continuous 
interpolation for p-degree basis functions with knot multiplicity k. Moreover, the ease in building spaces 
with high inter-element regularity allows for rather small problem sizes (in terms of degrees of freedom) 
with respect to standard finite element methods with the same approximation properties. This implies 

S^ | that, in general, for same approximation properties IGA stiffness and mass matrices are smaller than 

the corresponding finite element ones. Nevertheless, IGA discrete problems may still be very large in 
realistic problems of interest, and their condition numbers grow quickly with the inverse of mesh size h 
and the polynomial degree p. A detailed study of condition number estimates for stiffness matrix and 
mass matrix arising in isogeometric discretizations is given in [25]. For the /i-refinement, the condition 
number of the stiffness matrix is bounded above and below by a constant times h~ 2 , and the condition 
number of the mass matrix is uniformly bounded. For the p-refinement, the condition number is bounded 
above by p 2d 4P d and p 2 ^ d ~ 1 ^A pd for the stiffness matrix and the mass matrix, respectively. As a conse- 
quence, the cost of solving the linear system of equations arising from the isogeometric discretization, 
particularly using iterative solvers, becomes an important issue. Therefore, there is currently a growing 
interest in the design of efficient preconditioners for IGA discrete problems, in both the mathematical 
and the engineering communities. Multigrid methods for IGA have been introduced for two and three 
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dimensional elliptic problems by the authors in [24], and tearing and interconnecting methods for isoge- 
ometric analysis is discussed in [29]. Other recent work on solvers for IGA studied overlapping additive 
Schwarz methods [11, 13] and balancing domain decomposition by constraints methods [12]. Some 
issues arising in using direct solvers have been investigated in [18]. The results, we presented in [24], 
show optimal convergence rate with respect to the mesh parameter h. However, for discretizations based 
on higher degree polynomials, the convergence rate are quickly deteriorated. In this paper we discuss 
the construction of linear solvers which provide not only /i-independent convergence rates but also ex- 
hibit (almost) independence on p. The presented optimal order solvers are based on algebraic multilevel 
iteration (AMLI) methods. 

AMLI methods were introduced by Axelsson and Vassilevski in a series of papers [3, 4, 5, 6]. The 
AMLI methods, which are recursive extensions of two-level multigrid methods for FEM [2], have been 
extensively analyzed in the context of conforming and nonconforming FEM (including discontinuous 
Galerkin methods). For a detailed systematic exposition of AMLI methods, see the monograph [31, 38]. 
To reduce the overall complexity of AMLI methods (to achieve optimal computational complexity), var- 
ious stabilization techniques can be used. In the original work [3, 4], the stabilization was achieved by 
employing properly shifted and scaled Chebyshev polynomials. This approach requires the computation 
of polynomial coefficients which depends on the bounds of the eigenvalues of the preconditioned system. 
Alternatively, some inner iterations at coarse levels can be used to stabilize the outer iterations, which 
lead to parameter-free AMLI methods [5, 6, 30, 34]. These methods utilize a sequence of coarse-grid 
problems that are obtained from repeated application of a natural (and simple) hierarchical basis transfor- 
mation, which is computationally advantageous. Moreover, the underlying technique of these methods 
often requires only a few minor adjustments (mainly two-level hierarchical basis transformation) even if 
the underlying problem changes significantly. 

In this article we consider the scalar second order elliptic equation as our model problem. Let ft C 
M. d , d = 2,3, be an open, bounded and connected Lipschitz domain with Dirichlet boundary 9ft. We 
consider 

(1) - V • (AVu) = f in ft, u = u D on 9ft, 

where A{x) is a uniformly bounded function for x € ft. Let V° C -ff 1 (ft) denote the space of test 
functions which vanish on 9ft, and V D = V° + vP C H l (Q.) denote the set which contains the 
functions fulfilling the Dirichlet boundary condition on 9ft. By V® C V° and V® C V D we denote the 
finite-dimensional spaces of the B -spline (NURBS) basis functions. 
Introducing the bilinear form a(-, •) and the linear form /(•) as 



a(u,v) = / AWu-Vvdx, f(v)= / / 



(2) a(u, v) = I AVu-Vvdx, f(v)= / f v dx, 



the Galerkin formulation of this problem reads: 
Find Uh € V h D such that 

(3) a(u h , v h ) = f(v h ) for all v h £ Vfi. 

It is well known that (3) is a well-posed problem and has a unique solution. By approximating Uh and 
Vh using B -splines (NURBS) the variational formulation (3) is transformed in to a set of linear algebraic 
equations 

(4) An = f , 

where A denotes the stiffness matrix obtained from the bilinear form a(-, •), i.e. 

A= (oij) = (a(Ni,Nj)), i,j = l,2,3,....,n h , 

u denotes the vector of unknown degrees of freedom (DOF), and f denotes the right hand side (RHS) 
vector from the known data of the problem. Clearly, A is a real symmetric positive definite matrix. 

The rest of the paper is organized as follows. In Section 2 we briefly review the basics of B-splines 
and NURBS. An explicit representation of basis functions is also given in this section. The description 
of multilevel representation of B-splines (NURBS) is given in Section 3. A brief description of AMLI 
methods is given in Section 4. We then construct the isogeometric hierarchical spaces in Section 5. 
Numerical study of space splitting techniques is discussed in Section 6. The results of several numerical 



experiments in two- and three-dimensions are presented in Section 7. Finally, some conclusions are 
drawn in Section 8. 

2. B -SPLINES AND NURBS 

2.1. B-splines. We first recall the definition of B -splines, see e.g. [35,36]. 

Definition 1. Let Hi = {£j : i = 1, ..., n + p + 1} be a non-decreasing sequence of real numbers, 
called the knot vector, where £j is the i th knot, p is the polynomial degree, and n is the number of basis 
function. With a knot vector in hand, the B-spline basis functions, denoted by Nf(t;), are (recursively) 
defined starting with a piecewise constant 

(5a) 5i°(0 = { 1 m[6,6+l) * 

I otherwise, 

(5b) Bf (0 = /^-Sf- 1 ^) + / +P+1 ~/ £f+i(0, 

where < z < n,p > 1, owe? — jj considered as zero. 

The above expression is usually referred as the Cox-de Boor recursion formula, see e.g. [15]. For a 
B-spline basis function of degree p, an interior knot can be repeated at most p times, and the boundary 
knots can be repeated at most p + 1 times. A knot vector for which the two boundaiy knots are repeated 
p + 1 times is said to be open. In this case, the basis functions are interpolatory at the first and the last 
knot. Important properties of the B-spline basis functions include nonnegativity, partition of unity, local 
support and C p_fc -continuity. 

Definition 2. A B-spline curve C(£), is defined by 

(6) C7(£) = J>Bf (£) 

i=X 
where {Pi : i = 1, ..., n} are the control points and Bf are B-spline basis functions defined in (5). 

The previous definitions are easily generalized to the higher dimensional cases by means of tensor 
product. Using tensor product of one-dimensional B-spline functions, a B-spline surface £(£, rj) is de- 
fined as follows: 

(7) S^, V )=J2J2Bff^tr,)P i , j , 

t=l j=l 

where Pij,i = 1, 2, . . . ,m, j = 1, 2, . . . , ni, denote the control points, Bf 1 -' 1 ' 2 is the tensor product of B- 
spline basis functions Bf 1 and Bf, and Hi = {£i, £ 2 , . . . , £ ni+Pl+ i} and E 2 = {vi,m, • • • , Vn 2 + P2 +i} 
are the corresponding knot vectors. Similarly, B-spline solids can be defined by a three-dimensional 
tensor product. 

2.2. NURBS. While B-splines (polynomials) are flexible and have many nice properties for curve de- 
sign, they are also incapable of representing curves such as circles, ellipses, etc. exactly. Such limitations 
are overcome by NURBS functions. Rational representation of conies originates from projective geom- 
etry and requires additional parameters called weights, which we shall denote by w. Let {P™} be a set 
of control points for a projective B-spline curve in R 3 . For the desired NURBS curve in M 2 , the weights 
and the control points are derived by the relations 

(8) w i = (P* ) 3 , {P,i) d = {P?)/w t , d= 1,2, 

where Wi is called the i th weight and (Pi)d is the cP h -dimension component of the vector Pi. The weight 
function w(£) is defined as 



(9) ™(O = £a?(0 



i=i 



Then, the NURBS basis functions and curve are defined by 

The NURBS surfaces are analogously defined as follows 

m n 2 
i=l j=l 

where Nf 1 ^ 2 is the tensor product of NURBS basis functions Nf 1 and Nj 2 . Similarly, NURBS solids 
can be defined by a three-dimensional tensor product. NURBS functions also satisfy the properties of 
B-spline functions. For a detailed exposition see, e.g. [35, 36, 37]. 



2.3. Explicit Representation for B-splines. The recursive form of B-spline basis functions, given by 
(5), is elegant and concise, and is presented in all the IGA related references, see e.g., [15, 35, 36, 
37]. However, this form may not be the most efficient from computational point of view, specially 
when dealing with large knot vectors. To the best of authors' knowledge, within the IGA literature 
there is no reference on the explicit representation of B-splines for a given mesh size h. Therefore, 
we present the explicit form of B-splines in terms of the mesh size (knot-span) h. Having an explicit 
form of basis functions is also advantageous in devising inter-grid transfer operators for multigrid and 
multilevel iterative solvers. For brevity reasons, we restrict ourselves to a unit interval with equal spacing. 
Moreover, as most of the NURBS based designs in engineering use polynomial degree p = 2 and 3, we 
will confine ourselves up to p = 4 with C° and C p ~ 1 continuous basis functions. 



2.3.1. C p ~ 1 -continuity. We first consider the C p_1 continuous case as this is the default case for knot 
vector with non-repeated internal knots. For B-spline functions with p = and p = 1, we have the 
same representation as for standard piecewise constant and linear finite element functions, respectively. 
Quadratic B-spline basis functions, however, differ from their FEA counterparts. They are each identical 
but shifted related to each other, whereas the shape of a quadratic finite element function depends on 
whether it corresponds to an internal node or an end node. This "homogeneous" pattern continues for 
the B-splines with higher-degrees. 

We are interested to give an explicit representation for uniform B-spline basis functions defined on a 
knot vector E^ at any given level l^, k = 1, 2, 3, ..., with spacing h (= l/n), where n is the total number 
of knot spans. We shall use the notation Bf'[ for B-splines, where superscripts represent the polynomial 
degree and the regularity of basis functions, respectively, and the subscripts represent the level and the 
number of basis function, respectively. We start with level l\ with only one element. Using the definition 
from (5), at level l\ the B-spline basis functions of degree p = 2 on the knot vector E\ = {0, 0, 0, 1, 1, 1} 
are defined as follows 



(12) 



The mesh refinement takes place by inserting the knots. We consider uniform refinement of E\, i.e. 
inserting knots at the mid point of the knot values. At the next level I2, the basis functions on refined 
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knot vector Ei = {0, 0, 0, h, 1, 1, 1} are given by 



(13) 



B 



B 



B 



B 



2,p-l 

h,l 

2,p-l 
h,2 

2,p-l 
^2,3 

2,p-l 



'(l-2x) 2 , 0<x<±, 
0, \ < x < 1, 

'2x(2-3x), 0<x<±, 
2(1 -x) 2 , ~<x<l, 

'2x 2 , 0<x<±, 
-2 + 8x-6x 2 , ±<x<l, 

'0, 0<x<±, 
(l-2x) 2 , ±<x<l. 



Further refinements take place in a similar way, i.e., starting with E\, a single knot span, in the knot span 
Ek we will thus have 2 k ~ l knot spans. The explicit representation of B-splines at level If,, k > 3, is given 
by 



(14) 



B^^-ih-xf, 0<x<h, 



B 2,p-1 



B 



2,p-l 

Zfc,3+i 



— r^xUh — 3x), < x < h, 

2h z 

— -(2/i-x) 2 , h<x<2h, 

v 2/r 

^(x-zft) 2 , 

— + -(x - ih) - j^(x - ih) 2 , 



for/i< 



^(3h-(x-ih)) 2 



ih < x < (i + l)/i, 

(i + l)h<x< (i + 2)h, 

(i + 2)h <x < (i + 3)h, 
where i = 0, 1, 2, 3, ..., (l/7i) - 3, and /i < 1/4. 



B 2,p-1 
l k ,n+p-l 



— -2 (-1 + 2/1 + x) 2 , l-2/i<x < l-/i, 



-1 

12/? 



(3-4/i + 2(2/i-3)x + 3x 2 ), 1 - h < x < 1, 



for ft, < 



2' 



1 



^U^a-^-*) 2 , i-^<x<i. 



For higher degree polynomials, we can define the explicit representation in a similar way. Again using 
the definition (5) of B-splines, for p = 3, at first level l\ the basis functions with C p ~ l -continuity are 
given as follows 



B 



3,p-l 
Ji,l 



1-x) 



(15) 



S, 3 *" 1 = 3x(l - x) 2 

3,P-1 _ ^,3 



B 



h,2 

3,p-l 
!i,3 



3x 2 (l 



< x < 1, 
, < x < 1, 
, < x < 1, 



B 



hA 



0<x< 1, 



and at level I2, we have the following basis functions 



3, P -i_/(l-2x) 3 , 0<x<±, 



10, i <X< 1, 



I? 



■3>p-i _ (2x(3-9x + 7x 2 ), 0<x<i, 
'2,2 1 o/i _a3 1 «-- „, ^ 1 



2(1 -x) 3 , |<x<l, 



(16) 



3lP -i_/2x 2 (3-4x), 0<x<±, 
^2,3 



d3,P-1 
^2,4 



2(-l + x) 2 (-l + 4x), ±<x<l, 

2x 3 , 0<x<±, 

2- 12x + 24x 2 - 14x 3 , \ < x < 1, 



5; 



,3 1 p-l_JO, 0<x<± 



V 2 ,5 



2' 



-l + 2x) 3 , i < x < 1. 
For all other levels h,k > 3, the basis functions are denned below 



(17) 



Bf*' 1 = —(h-xf, 0<x<h, 



Bf' p ~ 1 = < 

'lb, 2 



x / 9x 7x 2 

1 



fix 2 



x 






6^ V Tft 



< x < /i, 
h < x <2h, 
0<x <h 



for /i < - . 

- 2' 



-3 9 



9 _ o 

T + Jjjr^ + Itf' ^^^ 

1 



777a (" 3/i + 



2h<x < 3/i, 



for /i < - , 



B 3,p-1 _ 



g3,p-l 



U/l 3 

r 1 o 

—r(x-ihY, ih<x<(i + l)h, 

bn 6 

2 2 1 1 

- - -(x - ih) + — 2 (x - i/j) 2 - — ^ (x - ihf, {i + l)<x<(i + 2)h, 

-22 10 4 1 

-7— + — (x - i/i) - -r(x - i/1) 2 + — t(x - i/1) 3 , (t + 2)h < x < (i + 3)/i, 
6 n h z Zh A 

f (l-^jr^) . (i + 3)fc<x<(i + 4)fc, 

where i = 0, 1, 2, 3, ..., (l//i) — 4, and h< -, 

r 7^3(-3^+(l-^)) 3 , 
-3 9(l-x) (1-x) 2 7(l-x) 3 

l (l-x) 2 / ll (l-x) 
,6 h? V 2 /i 



1 - 3h < x < 1 - 2h, 



l-h<x<l, 



f or /i < - 

~ 4 



d3,p-1 

Z fc ,n+p-l 



f ^~ 2h+il - X)) ^ 



— -^y — ^iL -r ys. — u.jj , 1 — 2h < x < 1 — h, 



/? 



/? 



/i 2 



f or /i < - , 
~ 2' 



B l P nl P = ^(h-{l-x))\ l-h<x<l. 



Finally, we give the explicit representation of basis functions for p = 4 with C p l -continuity. At level 
l\, with knot E\, the B-spline basis functions of degree p = 4 are given by 



(1,1 

ll,2 



{i- x y 



B7: v ~' = 4x(l - x) 



< x < 1, 



< x < 1, 



(18) 



B 
S 



4,p-l 
Zi,3 

4,p-l 



6x 2 (l-x) 2 , < x < 1, 



-4x 3 (l-x), 0<x<l, 
< x < 1. 



o4,p-l _ 4 
/i,5 ~~ X 



The B -splines on second level l<i with knot £?2 are defined as follows 



(19) 



B 



B 



B 



4,p-l 

4,p-l 
«2,2 

4,p-l 
«2,3 



r4, P -1 



o4,p-l 

^2,5 



o4,p-l 

^2,6 



;i - 2x) 4 , 

i 



0<x < \, 



0, i <X < 1, 

2x(4-18x + 28x 2 -15x 3 ), < x < \, 
2(1 -x) 4 , \<x<l, 

2x 2 (6-16x + llx 2 ), < x < \, 
2(1 -x) 3 (-l + 5x), l<x<l, 

2x 3 (4-5x), 0<x<±, 

2(1 -x) 2 (l-6x + llx 2 ), \<x<\, 

2x 4 , < x < §, 

-2 + 16a; - 48x 2 + 64x 3 - 30x 4 

0, 0<x<|, 

(l-2x) 4 , \<x<l. 



\ <x<l, 



At all other levels lk,k > 3, the basis functions of degree p = 4 with C p ^continuity are given by 



K7 1 = h* {h ~ x) *> °^ x<h > 



B 



4,p-l 

^,2 



-4a; 



/? 
1 



9x 7 x 2 15 x 3 
+ 4h~4h 2 + 32h 3 



(2h - x) 4 



B 



4,p-l 
*fc,3 



1 x 2 ( „ x 85 x 2 

9tf 27 - 33 S + W 



-3 x x 2 7x 3 23 x 4 



1 



U8/i 4 
( 2x 3 



5 



4,p-l 
Z fc ,4 



(3/i-x) 4 , 
25 x 4 



3 /i 3 72 /i 4 ' 

2 8x x 2 

3~3T + ¥ 



x 3 23 x 4 
2 h? + 72W 



22 40 x x 2 x 3 13 x 4 
r + Yh~ 8 h 2+2 h 3 ~72¥' 



24/i 4 



(4/i - x) 4 , 



< x < h, 
h < x < 2h, 

< x < h, 

h < x <2h, 
2h < x < 3/t, 

< x < /i, 
h< x <2h, 
2h < x < 3h, 
3h<x< Ah, 



for/i< -, 
~ 2' 



for h < -, 

~ 4 



for/i< 



(20) 






24/i 4 



(x - ih) 4 



ih < x < [i + l)h, 



1 I r 20 / -1.x 30 ^ -,n2 20 / -,n3 4 / -.N4 

— I -5 + — (x - ih) --fp{x- thy + -^ (x - thy - -^(x - thy 

(i + l)h < x < (i + 2)h, 

iA-2h (x - th) + m ix - th) -w {x ~ lh) -m {x - lh) > 

(i + 2)h < x < (i + 3)h, 

OOO 00 . ._ -. OO , -7 \2 " f -t \S / -7 \4 

— + -( x - lh )- w {x- l h) + ^{x-ih) - w (x-xh), 



1 



l24/i 4 



( 1 



B 



4,p-l 
Z fe ,n+p-3 



(5/i - (x - i/i)) - 



^-(1-x)) 4 



(i + 3)/i < x < (i + 4)h, 
(i + 4)h < x < (i + 5)h, 

where i = 0, 1, 2, 3, ..., (l//i) — 5, and h < -. 



2Ah A 

-22 40(1 -a;) Q (l-x) 2 



+ 



3 h h 2 

2 



+ 2 



1 - ih < x < 1 - 3/t, 

;i-x) 3 i3(i-x) 4 



13 



4,p-l 

fc ,n+p-2 



2 8(l-x) , (1-x) 

3 3 /i /i 2 

2 (l-x) 3 25 (1 - xf 
3 P 72 /? : 



( 3/l _ (1-x)) 4 , 



h 3 72 /i 4 ' 

1 - 3/i < x < 1 - 2/i, 

(l-x) 3 23 (l-x) 4 

/j3 + 72 /i 4 ' 

l-2/i<x<l-/i, 

l-/i <x < 1, 



18/i 4 



for h < - , 

1 - 3h < x < 1 - 2h, 



-3 , fi (l-s) 



(1-x) 2 7 (l-x) 3 23 (l-x) 4 
^ + 3 P 72 /i 4 : 



1 - 2/i < x < 1 - /i, 

1(1 -x) 2 / (1-x) 85(l-x) 2 \ 

' 27 - 33^ — '- + — V ,„ 7 ) , 1 - /i < x < 1, 



9 /i 2 



/? 



/ ? 2 



for h < - , 



B 4,p-1 

Z fc ,n+p-l 



(2/i-(l-x)) z 



1 

-4(1 -x 
h 



l-2h<x<l-h, 



9 (1 -a) 7 (l-x) 2 15 (l-x) 3 
+ 4 ft 4 /I 2 + 32 J? 

1 -h < x < 1, 



for/i< 



£ 



i,p-i 



1 



'fe,™+P /;4 



(/i-(l-x)) 4 , 0<x<l. 



2.3.2. C° -continuity. To lower the continuity of the basis functions across element boundaries, the knot 
values are repeated upto a desired level. By repeating the internal knots r times we get the C p ~ r continu- 
ous basis functions. In the previous section we have given the explicit representation for C p_1 continuous 
B-splines, which is the highest continuity for polynomial degree p. We now consider another extreme 
case, the lowest continuity, i.e. C° continuous basis functions. At first level l\ the C° continuous B- 
spline basis functions of degree p = 2, 3, 4 on a knot E\ = {0, 0, 0, 1, 1, 1} are same as those of C v ~ x 
continuous B -spline basis functions of same degree. 



The explicit representation for C° continuous B -spline basis functions of degree p = 2 at level l^,k > 
2 is given by 



(21) 



B 



2.0 
h 



ti = T2^- x >^ °< x <h, 



9 

B t ' 2+2i = -rfix — ih)(h + (x — ih)), (i — \)h < x < ih, 



B 



2.0 

i fe ,3+2i " 1 



p(/i + 0-i/i)) 2 



where i = 0, 1,2,3, ..., (l//i) - 1, 
(i — l)h < x < ih, 



(-h + (x-ih)) 2 , ih<x<(i + l)h, 

where i = 0, 1,2,3,..., {{1/h) -2) 



,2,0 



5^ = ^(1-/*-^, !-/»<*<!. 



For p = 3, the explicit representation for B-spline basis functions with C° -continuity, at level l^, k > 2, 
are given by 



(22) 



,3,0 



1 



B fc!i = p(>»-a0 3 > 0<*<fc, 



S 



3,0 



l k ,2+3i fo 



1 V 

1 + t(x — ih) I , ih < x < (i + l)h, 

where z = 0, 1, 2, 3, ..., (1/h) - 1, 



^ /-_ _-i.\2 ' - 



■Sj fc '" 3+ 3j = T2"(^ - i/i) I 1 - -(x - i/i) ) , ih < x < (i + l)h, 



^3,0 



ft 3 



(x - i/i) s 



where z = 0, 1, 2, 3, ..., (1/h) - 1, 
ih < x < (i + l)h, 



I -—(x-ih)) , (i + l)/t < x < (i + 2)h, 



where i = 0, 1, 2, 3, ..., ((l//i) - 2) , 



,3,0 = J_ 

l k ,np+l ^3' 



^(l-Zl-x) 3 , l-/l<X<l. 



Finally, the explicit representation for C° continuous basis functions of degree p = 4 at level If., k > 2, 
is given below 



B UA = hh-x)\ 0<x<h, 



(fc,l /;4 



B 



4,0 

Z fc ,2+4i — ^ 



-(x-tfe) 1 



(x — ih) 
h 



ih < x < (i + l)/i, 



where i = 0, 1, 2, 3, ..., (l//t) - 1, 



<'°3 + 4i = ^(^-^) 2 (l-^^) 



£ 



1,0 



,4+4* ^3 



where i = 0, 1, 2, 3, ..., (l/7i) - 1, 
(x-ihf (l- ( x ~ lh A ; i/j < x < (i + 1)/j, 
where i = 0, 1,2,3,..., (l//i) -1, 



B 



4,0 
l k ,5+ii 



(23) 



j^(x-ih) 4 , 



16 1 



2/; 



(x — ih) 



ih < x < (i + l)h, 

(i + l)h < x < (i + 2)h, 



where i = 0, 1,2,3, ...,((l//i) -2) 



s: 



1.0 



> 



h 



l k ,np+l- ~ X ) > "/1<X<1. 

3. Multilevel Representation of B-splines and NURBS 



3.1. Multilevel B-splines. In this section, we study the multilevel structure of B-splines and NURBS 
spaces. This will be used in the construction of corresponding hierarchical spaces (i.e. splitting the fine 
space into coarse space and its hierarchical complement) in Section 5. For a two level setting, let B p ' r 
and B p f denote the B -spline spaces at coarse and fine level, respectively. Let {B^'T, i = 1, 2, ..., n c } and 



{B 



p,r 



1, 2, ...,n./} be the set of basis functions for coarse and fine space, respectively, i.e. 



B p ' r 



spanl^,^,^,...,^;^}, 



and 



B 



p.r 



qT1aI1 (RP' r RP.' - nP. r RP' r \ 



The following result expresses coarse basis functions as the linear combination of fine basis functions. 



Proposition 3. Each coarse basis function B t 



p,r 



1, 2, ..., n c , can be represented as the linear com- 



p,r 



bination of the fine basis functions {B fi 
(24) 



1, 2, ..., rt/} by the following relation 



B v / 



oy By, 



i.e., 



B p ' r 

C.l 



X>^/j> 



where G p f = (gij)n c xn f , is called the restriction operator from a given fine level to the next coarse level 
for B -spline basis functions. 

In the following we explain the formation of transfer operator G p f at different levels of mesh and with 
increasing polynomial degree with both the extreme cases of C p ~ 1 and C° -continuity. 



3.1.1. C p ^-continuity. The B -spline basis functions B, 4 



2,p-l 



1,2,3, and B\ 



2,p-l 



->■! 



1,2,3,4, of 
degree p = 2 on knots E\ -- {0, 0, 0, 1, 1, 1} and E2 = {0, 0, 0, |, 1, 1, 1}, respectively, are defined in 
section 2.3.1. Clearly, the total number of coarse and fine basis functions are three (n c = 3) and four 
(rif = 4), respectively. The matrix G h ' p ~ = (gij)zx4 is given by the following representation of coarse 
basis functions as the linear combination of fine basis functions. 



D 2,p-1 o2,p-l . D 2.p-1 , D 2, 

B h,i = 5ii^ 2 'fi +5i2^/ 2 +giaB^ 



BZ' 1 = S^fS 1 + 9^B^ + g 23 B?^ + 9M Bf^\ 
Bf;^ 1 = aulffi- 1 + 932 BfS l + XxBfc 1 + 9uJB^\ 



Equivalently, it can be written as 
where 



B 



B 



2,p-l 
h 



B IT 

2,p-l 



B 



h,3 



r<2,p-i 



2,p-l _ ^.p-l^p-l 



511 512 913 914 
921 922 923 924 
931 932 933 934 



iBfr > 



bis 1 

<7 l 



For the above set of basis functions, G,' p is given by 



(25a) 



G?r 



4 2 
2 2 
2 4 



10 



Similarly, the coarse basis functions for B,' p - , i = 1, 2, 3, 4, at level li, can be obtained in terms of 



B 



2.P-1 



1, 2, ..., 6, by the following matrix 



(25b) 



Gir 



4 2 

2 3 10 

13 2 

2 4 



In a multilevel setting, the representation of each basis function B^ at level l^,k > 3, as the linear 

lk + 1 



combination of the basis functions B, ,p ~ . at level l^ + i is given by the the following matrix G h 



(25c) 



G 



2.P-1 



4 2 

2 3 1 

13 3 1 

1 3 3 



The size of the matrix G, ,p is 



3 3 1 

13 3 1 

1 3 2 

2 4 



ni k + 2) x {ni k+1 + 2), where ni k and raj fe+1 are the number of total 
l, respectively. 



knot spans at level If, and l^ 

For higher degree polynomials, the transfer operators can be defined in a similar way. For p = 3, 
at level l\ the basis functions B, ' p ~ ,i = 1,2,3,4, with C p ~ l -continuity can be represented by the 
following restriction operator 



(26a) 



The transfer operator for level 1$ can be written as 









2 


1 











G 3, P -1 = 


a) 


X 






1 



1 
1 




1 





















1 


2 



(26b) 



G 3, P -1 



For all levels l^+i with k > 3, we have 

16 8 



(26c) G 



3,p-l 

fc+l 



12 

4 



16 


8 




















8 


12 


3 

















4 


10 


4 

















3 


12 


8 




















8 


16 



3 
11 
2 



2 

12 



8 2 



The size of the matrix G z ' p is (m h + 3) x (nj fc+1 + 3). 



2 8 



12 


8 


2 




2 


8 


11 


4 






3 


12 



8 16 



ii 



Finally, we give the transfer operators for p = 4 with C p l -continuity. For levels I2 and Is the transfer 
operators are defined as follows: 



(27a) 









2 1 














/ 


i)* 


1 


1 











otr = ( 






1 




1 

1 1 


J 
















1 2 










' 48 24 















24 36 


9 








4,p-l 
h 


-U)« 


12 



30 

9 


9 
30 12 
















9 36 24 


















24 


48 



(27b) 



respectively. For levels l^, k > 3, the transfer operator is given by the following 

(27c) 

48 24 

24 36 9 

12 33 20 4 

6 25 29 15 3 

3 15 30 30 15 3 



g a, p - 



15 



30 
3 



30 
15 



15 

29 
4 



3 

25 
20 



6 
33 

9 



12 
36 



24 
24 



48 



where the size of the matrix is (m k + 4) x (ni k+1 + 4). 



Remark 4. Since in the span of an internal basis function of degree p at coarse level, there are p + 2 
full basis functions in the same span at fine level, therefore, any row of G^' p ~ can have at most p + 2 
entries. 



'fc+i 



3.1.2. C® -continuity. In section 2.3.2, we explained the explicit representation of C° continuous B- 
spline basis functions. The corresponding transfer operators are given in this section. The transfer 
2 with C°-continuity at level I2 is given by 



operator G, ' for p 



(28a) 



G 



2,0 



-,2,0 



The operator G l ,u for k > 2, in general, is given by 



(28b) G 



2,0 

fe+i 



2 1 
2 2 2 
1 2 



4 2 10 
2 2 2 








12 


4 


2 10 








2 2 2 
12 4 



12 



with size (2n; fc + 1) x (2rii h+1 + 1). The matrix G t ' , k > 2, has block structure with blocks G L ' 
The blocks are connected in such a way that if a block ends at ith row and jth column of G 



^,2,0 



2.0 
fc+i 



2 

then 



the next block will start at (i,j)th position of Gp with an overlap of last entry and first entry of the 
corresponding blocks. Note that, the first entry and the last entry in a block are same. 
The transfer operators for p = 3 with C°-continuity for levels li and 1% are given by 



(29a) 



and 



(29b) 



G 



3,0 



8 4 2 10 

4 4 3 2 

2 3 4 4 

12 4 8 



G 



3.0 



8 4 2 1 







4 4 3 


2 




2 3 
1 


4 4 







2 4 


8 


4 2 10 









4 4 3 2 






2 3 4 4 






12 4 8 



respectively. Following the same block structure as in G ; ' , we can generate G L ' with size (3n; fe + 
1) x (3n; fe+1 + 1). Finally for p = 4, we have the following transfer operators for levels I2 and I3 



(30a) 



and 



atf 



16 8421000 

8864200 

0466640 

0024688 

0001248 16 



(30b) 



G 



1.0 



16 8421000 




8864200 




0466640 




0024688 







12 4 8 


16 


8 4 2 10 







8864200 




0466640 




0024688 




0001248 16 



respectively. Similarly, repeating these blocks as in previous cases, we can generate G L ' with size 
(4n lh + 1) x (4n, fe+1 + 1). 

Remark 5. Note that the transfer operators are defined for one dimensional B-splines. For two- and 
three-dimensions, we take tensor product of these operators. 

3.2. Multilevel NURBS. This section presents the procedure for constructing NURBS multilevel spaces 
in a simplified manner. Since NURBS are generated from B-splines, its natural to construct NURBS 
transfer operators from B-splines transfer operators. For a two level setting, let A/J' r and A/?' r denote 



the NURBS spaces at coarse and fine level, respectively. Let {N t 



■p.r 



l,2,...,n c }and{iV 



1, 2, ...,n/} be the set of basis functions for coarse and fine space, respectively, i.e. 



p,r 

f,i> 



and 



Mr 



a/-- 



span{N%,N%,N%,...,N^ c }, 



span{N p f '[,N%,N : 



p.r 



N p ' r T 



13 



Note that, Proposition 3 holds for NURBS basis functions also, i.e., we have 



(31) M^ = RfMf r \ i.e., JV*;T = 5>tfJV£J, V* = l,2,3,..,n c> 

3=1 

where i?^' r = (r\y) ncX n/> is restriction operator with respect to NURBS basis functions. As NURBS are 
formed from B-splines and 'weights', RK' r can be obtained from G p / and 'weights'. Using the definition 
of NURBS and (31), we have 



?B p > r 



«^J 



(32) 



\ T ' A 



Vi = l,2,3,...,n c , 



i'=l 



i'=i 



where wf, i = 1,2, 3, ..., n c , and tc , j = 1, 2, 3, ...,n/, are the weights for coarse space and fine space, 

n 

respectively. Note that the weight function \, WiBi does not change its value with respect to refine- 
ments, i.e., we have 



i=l 



(33) 



E<^ = E^ p ' r 



f,r 



i=l j=l 

which is an important result from the refinement point of view. Now using (33), from (32) we get 

n f 



w 



3=1 



and thus 
(34) 



*sr = E 



n r f 



3=1 



W 



c /,3 



Comparing the coefficients of B v /- in (24) and (34), we get 



(35) 

This can be equivalently written as 
(36) 



r ij w j 



9ij 



r _ 9iJ w i 
<ij — * ■ 



TC n P,r f-fTrf 



Ry =wfGy(w{ 



where Wj and Wj are the diagonal matrices corresponding to weights at the coarse level and fine level, 
respectively, and defined as follows 



Wf 



w\ 



w:-, 



w 



n c —l 



w; 



,wj 



UK 



(Co 



w 



nt—1 



W 



f 



The equation (36) gives us the relation between B-splines and NURBS operators using B-splines transfer 
operators and weights at coarse and fine levels. From (33) we can also obtain the procedure to refine the 
weights as follows. We have 



"/ 



E< B c;[ = E^ r 



fj> 



i=l 



3=1 
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which implies 

n c n s n f 

Comparing the coefficients of B v f- from both the sides, we get 



(37) 



W J 



j =^2 w i9ij for j = 1,2, ..., 



rn 



i=\ 



Equivalently, the above can be written in matrix form as follows 



(38) 
where 



W f == (Gf 



W c 



W f 



UK 



if:-., 



IV 



nr — 1 



IV 



f 



H J 



,w c 



iiA 



w 



n c —l 



Using above, now we can write the NURBS operators in terms of B -spline operator and weights only at 
coarse level. From (35), we get 

(39) 



i=l 

In matrix form the above can be written as 

(40) Rf = WfGf (diag ( (Gf ) ' W c 



T 



Remark 6. The operators G p / and R p / can also be used in constructing restriction operators in multi- 
grid methods, see e.g., [24]. 

4. AMLI Methods 

In this section we present the basic principle of AMLI methods. In what follows we will denote by 
M^ fc ) a preconditioner for stiffness matrix A^ lk ' corresponding to level If.. We will also make use of the 
corresponding hierarchical matrix A^ lk \ which is related to A^ lk > via a two-level hierarchical basis (HB) 
transformation J^ k > , i.e., 
(41) i(W = J (W,4fe)(jfe))T 

The transformation matrix J^ lk ' specifies the space splitting, which will be described in detail in Sec- 
tion 5. By A]- and A] , 1 < i,j < 2, we denote the blocks of A^ and A^ that correspond to the 
fine-coarse partitioning of degrees of freedom (DOF) where the DOF associated with the coarse mesh 
are numbered last. 

The aim is to build a multilevel preconditioner M^*> for the coefficient matrix A^> := Ah at the level 
of the finest mesh that has a uniformly bounded (relative) condition number 

x(mW _1 aW) = 0(1), 

and an optimal computational complexity, that is, linear in the number of degrees of freedom Nf at the 
finest mesh (grid). In order to achieve this goal hierarchical basis methods can be combined with various 
types of stabilization techniques. 

One particular purely algebraic stabilization technique is the so-called algebraic multilevel iteration 
(AMLI) method, where a specially constructed matrix polynomial p^ of degree v^ can be employed at 
some (or all) levels l^. The AMLI algorithm has been originally introduced and studied in a multiplicative 
form, see [3, 4]. 
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We have the following two-level hierarchical basis representation at level l^ 



(42) 



AM 



2(l k ) 



A 



(W 

12 



A 



(h) j(i») 



21 



.4 



22 



.4 



i(lk) 



' A ih) 

4i } ^ ( ' fe - i} 



12 



Starting at level l\ (associated with the coarsest mesh), on which a complete LU factorization of the 
matrix A^ ll > is performed, we define 

(43) M (/l) := A {h) . 

Given the preconditioner M( lk ~ 1 ' at level h-i, the preconditioner M" k > at level Ik is then defined by 

(44) M {h) .- L {i k )ij{h) ^ 

where 



(45) 



l (!h) ... 



r (h) 

pk) 



a 



o 
(h) 

22 



£/(**) .- 



c\ 



(hT l 2(J fc ) 



.4 



12 



Here C}^ is a preconditioner for the pivot block -A^ * , and 



(46) 



C 



C»fc) 

22 



(47) 0<p (w (t)<l, 

It is easily seen that (46) is equivalent to 



ph-i) (i-p^^M^-^ 1 A {lk - l) )\ 



< t < 1, p ( ' fc) (0) = l. 



(48) 



C 



22 



^ffe-l) 1 g('fc)(^4( / fc-l)_/\^(^-l) 1 ) 



where the polynomial q" k ' is given by 
(49) 



M) 



(x) 



1 -p^Xx) 



We note that the multilevel preconditioner defined via (44) is getting close to a two4evel method when 



q" k >{x) closely approximates 1/x, in which case C\^ 

y(h) 



(hY 



A (h 



. In order to construct an efficient 



multilevel method, the action of C^ 2 on an arbitrary vector should be much cheaper to compute (in 

terms of the number of arithmetic operations) than the action of A^- 1 ' . Optimal order solution algo- 

(l )~ l 
rithms typically require that the arithmetic work for one application of C^ is of the order 0(Ni k _ 1 ) 

where Ni k _ 1 denotes the number of unknowns at level lk-i. 

It is well known from the theory introduced in [3, 4] that a properly shifted and scaled Chebyshev 

polynomial p^ 1 ^ := p Vh of degree Vk can be used to stabilize the condition number of M^ k ' A^ lk ' (and 
thus obtain optimal order computational complexity). Other polynomials such as the best polynomial 
approximation of l/x in uniform norm also qualify for stabilization, see, e.g., [32]. Alternatively, in the 
nonlinear AMLI method, see, e.g., [6], a few inner flexible conjugate gradient (FCG) type iterations (for 
the FCG algorithm, see also [33]) are performed in order to improve (or freeze) the residual reduction 
factor of the outer FCG iteration. In general, the resulting nonlinear (variable step) multilevel precon- 
ditioning method is almost equally efficient, and, because its realization does not rely on any spectral 
bounds, is easier to implement than the linear AMLI method (based on a stabilization polynomial). For 
a convergence analysis of nonlinear AMLI see, e.g., [30, 38]. 

The construction of optimal preconditioners in the framework of AMLI methods is based upon a 
theory in which the constant 7 in the strengthened Cauchy-Bunyakowski-Schwarz (CBS) inequality plays 
a key role. The CBS constant measures the cosine of the abstract angle between the coarse space and 
its hierarchical complementary space. The general idea is to construct a proper splitting by means of a 
hierarchical basis transformation. 

In the hierarchical bases context we denote by V\ and V% subspaces of the space V/ t . The space Vi is 
spanned by the coarse-space basis functions and V\ is the complement of V2 in Vh , i.e., Vh is a direct 
sum of V\ and Vi : 

V h = Vi © v 2 . 
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Let Vi € Vi, i = 1,2. The CBS constant measures the strength of the off-diagonal blocks in relation 
to the diagonal blocks (see, (42)) and can be denned as the minimal 7 satisfying the strengthened CBS 
inequality 



(50) 



\vfAi 2 v 2 \ < r yUvTAiiv 1 )(v 2 r A 2 2V2)\ 



1/2 



A detailed exposition about the role of this constant can be found in [22]. 

Typically, the iterative solution process is of optimal order of computational complexity if the degree 
u k = v of the matrix polynomial (or alternatively, the number of inner iterations for nonlinear AMLI) at 
level Z/c satisfies the optimality condition 



(51) 



l/V(l-7 2 )< »< r, 



where r ps 77. = Ni k /Ni k _ x denotes the reduction factor of the number of degrees of freedom (DOF), 
and 7 denotes the constant in the strengthened Cauchy-Bunyakowski-Schwarz (CBS) inequality. The 
value of r is approximately 4 and 8 in case of two- and three-dimensional problems, respectively. For a 
more detailed discussion of AMLI methods, including implementation issues, see, e.g., [31, 38]. 

Remark 7. The preconditioner defined in (44) is of multiplicative form. The introduction of AMLI meth- 
ods was based on the multiplicative form, see [3, 4, 5, 6], and is commonly used in practice. However, it 
is also possible to choose the preconditioner in the additive form, which is defined as follows 



(52) 



M 



(lk) .. 



a 





22 



In this case the optimal order of computational complexity demands that the matrix polynomial degree 
(or the number of inner iterations of nonlinear AMLI) satisfy the following relation 



(53) 



V(l + 7)/(l-7) < v < t. 



5. Construction of Hierarchical Spaces 

The hierarchical basis techniques result in splittings in which the angle between the coarse space and 
its hierarchical complement is uniformly bounded with respect to the mesh size. Recall from section 4, 
we have the following two-level hierarchical basis representation for stiffness matrix at fine level 



Aif) 






A 



if) 

21 



M)' 

m 

^22 



M) 



.4 



21 



A® 



"(f) "(f) 

where A 22 represents the matrix corresponding to coarse basis functions and A\\ represents the matrix 

corresponding to its hierarchical complement. Recall from Section 3, for B -splines and NURBS we have 

the following transformations 

l(/) _ a(c) rtf>,r Af(/~iP,r\T 
A (c) 



(54) 



A 
A 



22 



22 



GfAf(G p / j 
R p fAf(EPff, 



respectively. For hierarchical complementary spaces, let T p / be the matrix such that 



/ 



(55) 



A 



(/) 
11 



T P>r A f( T P,r^T 



7 



7 



Here the matrix T P ' r is a hierarchical complementary transfer operator, which transfers fine basis func- 
tions to a set of hierarchical complementary basis functions. The remaining two blocks of the hierarchical 
matrix AV' can be obtained by the following relations 

A^ =Tf r AW(G p ff, 

* U) - G p ; r A^{T P /) T . 



21 



To construct T P,r efficiently, the following points are important. 



/ 



m 



(1) The condition number of the AW block should be independent of mesh size. 
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(2) The CBS constant 7, see (50), should be bounded away from one, i.e. the minimum generalized 
eigenvalue of block A^ 2 with respect to the Schur complement should be greater than 1/4 for 
v = 2 and 1/9 for v = 3. 

(3) The basis for hierarchical complementary space should be locally supported. In other words, the 
block An should be sparse and nicely structured. 

The construction of T P,T , based on the linear combination of fine basis functions, is not unique. Based 
on the above mentioned guidelines, a representation of a complementary basis function should not in- 
volve several fine basis functions, because it will cause less sparse structure of matrix T P ' r . Based on 
our extensive study with different choices of linear combinations satisfying the above requirements, we 
present two choices for T P,r , for p = 2, 3, 4 and for extreme cases of smoothness, namely C p_1 and C . 
For the first choice, we have the following matrix representation of hierarchical complementary space 
imp = 2 with C p_1 -continuity. 



->2,p-l 
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-1 
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-1 










1 


-1 












1 


-1 



1 


-1 








1 


-1 



The above matrix has the block structure with blocks, say M 1 ' p . The blocks are connected in such a 
way that if a block ends at ith row and jth column of T, ' p ~ then the next block will start at (i+ 1, j — l)th 

position of T L ,p ~ . In general, for p = 2, 3, 4, we write the following block form of T p ' p ~ with blocks 
m p,p-i 



(57) 



where 



rpP,P~l 
lk+1 



M p 



p-i 



m p,p-i 



m p,p-i 



M{ 



p,p-i 



and 



M 


2,p-l = 


M^ 1 = 


"0 -1, 



M^ 1 = 


' 1/2 




1 




-1 

1 




-1 



■1/2 3/4 -1/2 

-1/2 3/4 -1/2 



-1 



1 -1/2 
1/2 -1 1 -1/2 

respectively. The blocks are connected in such a way that if a block ends at ith row and jth column of 
Tf ,p ~ l then the next block will start at (i + 1, j - (p - l))th position of Tf ,p ~ 1 . 
For second choice of T P,r we give the following block matrix 



(58) 





- M™' 1 




" 


rpp,p-l _ 
lk+1 


M™- 1 


M™- 1 






M™" 1 J 


'~ are given by 




M^- 1 


= 


' -1/2 1 -1 
-1/2 


1/2 
1 -1 1/2 


1 



IK 



and 



M 



M 



4,p-l 



3.P-1 



1/8 -1/2 3/4 -1/2 1/8 
1/8 -1/2 3/4 -1/2 1/8 



1/4 1/2 -1 1 -1/2 -1/4 
1/4 1/2 -1 1 -1/2 -1/4 



respectively. 

For C° continuous basis functions, we give the following matrix representation of hierarchical com- 
plementary spaces. For p = 2 we have 



^2,0 
lk+1 



1 -1/4 
1 -1/4 



1 -1/4 
1 -1/4 



1 




-1/4 
1 





-1/4 



The above matrix has the block structure and the blocks are connected in such a way that if a block ends 
at ith row and jth column of T x ' then the next block will start at (i + 1, j)th position of 7} ' . In 
general, for p = 2, 3, 4, we can write the following hierarchical complementary operators 



(59) 



where 



-,p,0 

lk + 1 



M- 



1>S) 



M- 



p,a 



Mf' 1 



Mf 



M 



2,0 



1 




-1/4 
-1/4 1 



M- 



3,0 



1 













1/2 






-1/2 
1 










1 



and 



M 



1,0 



-2/3 5/4 

-2/3 5/4 

5/4 








-2/3 
5/4 















-2/3 



respectively. Another choice for T?' r for C° continuous basis functions is obtained by choosing the 



following block matrix 



7 



(60) 



where 



T Pfi 

( fc+i 



Mr 



M. 



3,0 





" Mf' u 






Mf 






Mf 


1 




Mf 


- 


!,0 


' -1/4 1 -1/4 






-1/4 1 -1/4 


; 


-1/2 1/2 





-1/4 1/10 -1/4 








1/2 -1 


/2 
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and 



ur = 


" 





-5/9 





1 

-5/9 




-5/9 
1 






-5/9 

-5/9 







1 

-5/9 





-5/9 
1 






-5/9 


respectively. 



















Remark 8. All the above operators are defined for one-dimension. The higher dimensional operators 
are obtained via tensor product. 

6. Numerical results 

To test the performance of the AMLI methods, we consider the following test problems, whose dis- 
cretizations are performed using the Matlab toolbox GeoPDEs [21, 23]. 

Example 9. Let = (0, l) 2 . Together with A = I, and Dirichlet boundary conditions, the right hand 
side function f is chosen such that the analytical solution of the problem is given by u = e x sin(y). 

Example 10. The domain is chosen as a quarter annulus in the first Cartesian quadrant with inner 
radius 1 and outer radius 2. Together with A = I, and homogeneous Dirichlet boundary conditions, the 
right hand side function f is chosen such that the analytic solution is given by u = —xy 2 (x 2 + y 2 — 
l)0 2 + y 2 -4), see [21,23]. 

Example 11. The domain is chosen as a thick quarter of a ring. Together with A = I, and Dirichlet 
boundary conditions, the right hand side function f is chosen such that the analytical solution of the 
problem is given by u = e x sm(xy)cos(z). 

At the finest level (largest problem size), the parametric domain is divided into n equal elements in 
each direction. The initial guess for (iteratively) solving the linear system of equations is chosen as 
the zero vector. Let tq denote the initial residual vector and ru denote the residual vector at a given 
PCG/FCG iteration n; t . The following stopping criteria is used 

(61) Jp4 < !0~ 8 - 

Foil 

/ Until \V n it 
The average convergence factor reported in the following tables is defined as p = - — - . In 

the following tables, by LI, L2 and N2 we denote the linear multiplicative AMLI cycles with v = 1, 
v = 2 and non-linear multiplicative AMLI cycle with v = 2, respectively. By t c , we represent the 
setup time, i.e., the time taken in the construction of transfer operators and generating the preconditioner 
for An block (for which we used the ILU(O) factorization, i.e. without any fill-in). The solver time is 
represented by t s . For all the test cases we take the polynomial degree p = 2,3,4 with C°- and C p_1 - 
continuity. Furthermore, the transfer operator Gf' r is fixed and it exactly represents the coarse basis 

functions in the space of fine basis functions. The hierarchical complementary transfer operator Tf' r 
are chosen in two different ways as defined in Section 5, see (57)-(60). 

We first consider the Example 1 and provide t c , t s , mt and p for L1-, L2-, N2- cycles with both the 
choices of Tf' r . Numerical results are presented in Tables 1-2 and Tables 3-4 for first choice and second 

choice of Tf ,r , respectively . From Tables 1-4 we observe the following: 

• The number of iterations and total solution (t c + t s ) time show an /i-independent convergence 
rates for C p ~ 1 - and C°-continuity. 

• For C p_1 -continuity, the results are almost p-independent, whereas for C° -continuity, the degree 
p has some effect on PCG/FCG iterations. 

• For C p ~ l -continuity, all the AMLI cycles give optimal results, and the V-cycle (y = 1) is the 
fastest among all. This is due to a very nice bound on 7 for C p_1 -continuity. Therefore, in the 
remaining numerical computations we consider linear AMLI cycle with v = 1 and nonlinear 
AMLI cycle with v = 2 for C p ~ l continuous basis functions. 

• For C° -continuity, V-cycle [y = 1) is not an optimal order method, an observation similar to 
standard FEM. However, for C° -continuity, both the v = 2 cycle methods (linear and nonlinear) 
exhibit optimal order behavior, and nonlinear AMLI gives overall better results. Therefore, we 
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Table 1. AMLI methods for Example 1: First choice of Tf ,r given in (57) with C p 1 regularity 



1/h 


te 


t. 


n it 


P 






LI 


L2 


N2 


LI 


L2 


N2 


LI 


L2 


N2 












p = 


2 










8 


0.00 


0.00 


0.00 


0.00 


1 


7 


7 


0.0641 


0.0641 


0.0622 


16 


0.00 


0.00 


0.01 


0.01 


8 


7 


7 


0.0948 


0.0966 


0.0670 


32 


0.01 


0.01 


0.01 


0.01 


9 


8 


7 


0.1108 


0.0988 


0.0672 


64 


0.04 


0.02 


0.03 


0.04 


9 


8 


7 


0.1086 


0.0901 


0.0622 


128 


0.18 


0.07 


0.09 


0.12 


9 


8 


7 


0.1166 


0.0909 


0.0624 


256 


0.72 


0.25 


0.30 


0.41 


9 


8 


7 


0.1175 


0.0879 


0.0603 


512 


2.97 


1.03 


1.12 


1.50 


9 


8 


7 


0.1276 


0.0945 


0.0620 












P = 


3 










8 


0.00 


0.00 


0.00 


0.00 


8 


8 


8 


0.0901 


0.0901 


0.0901 


16 


0.01 


0.01 


0.01 


0.01 


9 


9 


8 


0.1111 


0.1129 


0.0686 


32 


0.02 


0.01 


0.01 


0.02 


10 


9 


7 


0.1293 


0.1043 


0.0577 


64 


0.10 


0.03 


0.04 


0.05 


10 


8 


7 


0.1361 


0.0857 


0.0551 


128 


0.41 


0.12 


0.13 


0.18 


10 


8 


7 


0.1369 


0.0821 


0.0536 


256 


1.76 


0.48 


0.46 


0.63 


10 


8 


7 


0.1348 


0.0794 


0.0523 


512 


7.50 


1.65 


1.77 


2.37 


9 


8 


7 


0.1283 


0.0771 


0.0511 












P = 


4 










8 


0.00 


0.00 


0.00 


0.00 


10 


10 


10 


0.1139 


0.1139 


0.1139 


16 


0.01 


0.01 


0.01 


0.01 


12 


12 


10 


0.1866 


0.1882 


0.1378 


32 


0.06 


0.02 


0.02 


0.03 


12 


11 


9 


0.2013 


0.1822 


0.1100 


64 


0.26 


0.07 


0.07 


0.10 


12 


10 


9 


0.2038 


0.1557 


0.1032 


128 


1.09 


0.26 


0.24 


0.37 


12 


9 


9 


0.2028 


0.1209 


0.0977 


256 


4.57 


0.98 


0.88 


1.21 


12 


9 


8 


0.1976 


0.1182 


0.0975 


512 


19.05 


3.60 


3.44 


4.66 


11 


9 


8 


0.1853 


0.1146 


0.0930 



consider only nonlinear AMLI cycle with v = 2 for C° continuous basis functions in remaining 

numerical results. 

For p = 4 with C p ~ 1 -continuity, we could not obtain better 7 with the second choice of Tf ,r as 

compared to the first choice. Therefore, in Table 3, the numerical results are presented only for 

p = 2, 3 with second choice of T?' r . Numerical results for p = 4 may be improved by choosing 

different operators, which demands further investigation. 

For C p_1 -continuity, though the number of iterations are less for second choice of Tf ,r , the 

overall time (t c +t s ) is more than the first choice of T^' r . This happens due to comparatively less 



sparse structure of second choice Tf ,T , which results in more construction time t c . Therefore, 
in the remaining numerical tests we consider only the first choice of Tf' r for C p_1 continuous 
basis functions. 

For C°-continuity, we get mixed results from both the choices of T[' r . This is due to the fact 
that there is not much difference in number of nonzero entries in each row of T?' r for two 



different choices. Therefore, numerical results are provided for both the choices of T? ,r 



forC 

continuous basis functions. 
We now consider Example 2 with curved boundaiy. The geometry for this example is represented 
by NURBS basis functions of order 1 in the radial direction and of order 2 in the angular direction, see 
[21]. Numerical results are provided for C p_1 -continuity with first choice of Tf ,r in Table 5, and for 

C°-continuity with both the choices of Tf' r in Table 6. All the results are qualitatively similar to that of 
Example 1 with square domain. 

Finally, we consider three-dimensional problem as stated in Example 3. The numerical results are pre- 
sented in Tables 7-8. Due to the limitation of computer resources numerical results in three-dimensions 



*did not converge upto desired accuracy. 
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Table 2. AMLI methods for Example 1: First choice of Tf ,r given in (59) with C° regularity 



1/h 


tc 


t s 


n it 


P 






LI 


L2 


N2 


LI 


L2 


N2 


LI 


L2 


N2 


p = 2 


8 


0.00 


0.01 


0.01 


0.01 


9 


9 


9 


0.1072 


0.1072 


0.1072 


16 


0.01 


0.01 


0.01 


0.01 


11 


11 


9 


0.1695 


0.1716 


0.1102 


32 


0.02 


0.03 


0.04 


0.04 


13 


11 


9 


02195 


0.1738 


0.1110 


64 


0.10 


0.07 


0.09 


0.10 


14 


11 


9 


0.2606 


0.1744 


0.1109 


128 


0.38 


0.30 


0.29 


0.34 


16 


11 


9 


0.2973 


0.1743 


0.1105 


256 


1.65 


1.25 


1.04 


1.23 


17 


11 


9 


0.3288 


0.1736 


0.1102 


512 


6.93 


5.17 


3.84 


4.61 


18 


11 


9 


0.3557 


0.1730 


0.1100 


p = 3 


8 


0.01 


0.01 


0.01 


0.03 


12 


12 


12 


0.1999 


0.1999 


0.1999 


16 


0.02 


0.03 


0.04 


0.03 


17 


17 


12 


0.3288 


0.3305 


0.2124 


32 


0.09 


0.09 


0.12 


0.10 


22 


18 


12 


0.4258 


0.3568 


0.2129 


64 


0.37 


0.38 


0.41 


0.33 


27 


19 


12 


0.5014 


0.3650 


0.2122 


128 


1.55 


1.77 


1.34 


1.21 


32 


19 


12 


0.5581 


0.3673 


0.2114 


256 


6.73 


7.87 


4.88 


4.51 


37 


19 


12 


0.6038 


0.3670 


0.2110 


512 


28.76 


36.56 


19.12 


17.84 


42 


19 


12 


0.6394 


0.3664 


0.2108 


p = 4 


8 


0.01 


0.03 


0.03 


0.03 


19 


19 


19 


0.3631 


0.3631 


0.3631 


16 


0.05 


0.07 


0.10 


0.09 


25 


26 


19 


0.4784 


0.4827 


0.3719 


32 


0.24 


0.32 


0.36 


0.29 


38 


30 


19 


0.6087 


0.5337 


0.3720 


64 


1.07 


1.62 


1.29 


1.05 


52 


32 


19 


0.6982 


0.5585 


0.3719 


128 


4.50 


8.04 


4.90 


3.89 


67 


34 


19 


0.7585 


0.5766 


0.3709 


256 


18.73 


40.11 


19.41 


15.25 


85 


35 


19 


0.8038 


0.5827 


0.3703 


512 


76.22 


190.29 


77.37 


62.24 


100* 


35 


19 


0.8379 


0.5878 


0.3700 



Table 3. AMLI methods for Example 1: Second choice of 7] p,r given in (58) with 
QP- 1 regularity 



1/h 


tc 


t a 


n it 


P 






LI 


L2 


N2 


LI 


L2 


N2 


LI 


L2 


N2 












p = 


2 










8 


0.08 


0.02 


0.42 


0.52 


5 


5 


5 


0.0227 


0.0227 


0.0227 


16 


0.00 


0.01 


0.01 


0.01 


6 


6 


5 


0.0304 


0.0326 


0.0217 


32 


0.02 


0.01 


0.01 


0.01 


6 


6 


5 


0.0316 


0.0311 


0.0226 


64 


0.07 


0.02 


0.05 


0.05 


6 


6 


5 


0.0303 


0.0300 


0.0224 


128 


0.30 


0.06 


0.08 


0.10 


6 


6 


5 


0.0314 


0.0310 


0.0234 


256 


1.21 


0.22 


0.30 


0.39 


6 


6 


5 


0.0301 


0.0296 


0.0226 


512 


5.18 


0.88 


1.05 


1.62 


6 


6 


6 


0.0326 


0.0321 


0.0269 












p = 


3 










8 


0.00 


0.02 


0.00 


0.00 


1 


7 


7 


0.0443 


0.0443 


0.0443 


16 


0.01 


0.00 


0.00 


0.01 


1 


7 


6 


0.0560 


0.0569 


0.0365 


32 


0.04 


0.01 


0.01 


0.02 


1 


7 


6 


0.0576 


0.0494 


0.0319 


64 


0.18 


0.04 


0.04 


0.05 


1 


6 


5 


0.0569 


0.0377 


0.0216 


128 


0.77 


0.12 


0.13 


0.17 


1 


6 


5 


0.0542 


0.0343 


0.0204 


256 


3.32 


0.47 


0.49 


0.64 


1 


6 


5 


0.0502 


0.0326 


0.0195 


512 


13.99 


1.60 


1.89 


2.44 


6 


6 


5 


0.0446 


0.0311 


0.0186 



are provided only upto h 
with v 



1/32. In Table 7, linear AMLI cycle with v = 1, and nonlinear AMLI cycle 



2 are given for C p 1 continuity with first choice of T?' 1 . The results exhibit optimal order for 
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Table 4. AMLI methods for Example 1: Second choice of Tf ,r given in (60) with C° regularity 



1/h 


*c 


ts 


nu 


P 






LI 


L2 


N2 


LI 


L2 


N2 


LI 


L2 


N2 










I 


= 2 












8 


0.00 


0.00 


0.00 


0.00 


8 


8 


8 


0.0901 


0.0901 


0.0901 


16 


0.01 


0.01 


0.01 


0.01 


9 


9 


8 


0.1173 


0.1195 


0.0918 


32 


0.04 


0.02 


0.03 


0.04 


10 


9 


8 


0.1308 


0.1197 


0.0900 


64 


0.15 


0.07 


0.09 


0.12 


10 


9 


8 


0.1430 


0.1192 


0.0890 


128 


0.62 


0.27 


0.32 


0.41 


10 


9 


8 


0.1479 


0.1191 


0.0884 


256 


2.71 


1.08 


1.20 


1.56 


10 


9 


8 


0.1509 


0.1191 


0.0880 


512 


11.22 


4.33 


4.55 


6.00 


10 


9 


8 


0.1528 


0.1191 


0.0878 










I 


= 3 












8 


0.01 


0.01 


0.03 


0.03 


9 


9 


9 


0.1133 


0.1133 


0.1133 


16 


0.02 


0.02 


0.02 


0.02 


11 


11 


9 


0.1724 


0.1743 


0.1191 


32 


0.09 


0.05 


0.07 


0.07 


13 


11 


9 


0.2216 


0.1762 


0.1206 


64 


0.39 


0.20 


0.22 


0.25 


14 


11 


9 


0.2627 


0.1777 


0.1212 


128 


1.62 


0.88 


0.79 


0.91 


16 


11 


9 


0.2998 


0.1782 


0.1215 


256 


7.06 


3.70 


2.87 


3.42 


17 


11 


9 


0.3321 


0.1785 


0.1216 


512 


29.78 


16.54 


11.24 


13.37 


19 


11 


9 


0.3630 


0.1786 


0.1217 










I 


= 4 












8 


0.01 


0.02 


0.02 


0.02 


10 


10 


10 


0.1368 


0.1368 


0.1368 


16 


0.07 


0.04 


0.06 


0.05 


13 


13 


10 


0.2199 


0.2219 


0.1419 


32 


0.33 


0.15 


0.18 


0.18 


15 


13 


10 


0.2825 


0.2254 


0.1416 


64 


1.88 


0.61 


0.60 


0.64 


17 


13 


10 


0.3259 


0.2251 


0.1412 


128 


5.92 


2.51 


2.20 


2.41 


18 


13 


10 


0.3575 


0.2247 


0.1410 


256 


25.51 


11.39 


8.56 


9.56 


20 


13 


10 


0.3844 


0.2245 


0.1408 


512 


104.33 


49.49 


35.15 


39.67 


21 


13 


10 


0.4042 


0.2244 


0.1408 



both the solvers. The increased number of iterations (as compared to two-dimensional examples) can be 
attributed to the smaller angle between coarse space and its complementary space. For C°-continuity the 
numerical results with both the choices of T? ,r are given in Table 8. The first choice of Tf '' , however, 
does not result in an optimal order method. The optimality is restored with v = 3, which are presented in 
the column with N3. The second choice, though expensive, gives optimal order method for second order 
stabilization [v = 2). In Tables 7-8, The entries marked by * represent the cases where the computations 
are performed on a machine with larger memory but shared with other users, therefore timings are not 
provided for these cases. 

We note that for two-dimensional problems, the total time of the solvers also exhibit optimal com- 
plexity, however, for three-dimensional problem the increase in the total time (t c + t s ) for successive 
refinement is more than the factor of increase in number of unknowns. This is due to the construction 
of operators Gf ,r and 1? ,r by tensor product of matrices for one-dimensional operators (see Remark 

8), and expensive preconditioner for An (ILU(0)). In our future study on local analysis, we also intend 
to construct these operators for two- and three-dimensional problems without tensor product, and devise 
efficient and cheaper preconditioner for A\\. 

1. Conclusions 

We have presented AMLI methods for the linear system arising from the isogeometric discretization of 
the scalar second order elliptic problems. We summarize the main contribution of this paper as follows. 

(1) We provide the explicit representation of B -splines as a function of mesh size h on a unit interval 
with uniform refinement. The explicit representation is given for C° and C p ~ l continuous basis 
functions of polynomial degree p = 2, 3, 4, the most widely used cases in engineering applica- 
tions. Explicit form of B-splines is important from computational point of view, as well as in 
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Table 5. AMLI methods for Example 2: First choice of 7] p ' r given in (57) with C p x regularity 



'fe+i 



1/h 


tc 


ts 


nu 


P 






LI 


N2 


LI 


N2 


LI 


N2 


p = 2 


8 


0.02 


0.02 


0.01 


8 


8 


0.0802 


0.0802 


16 


0.00 


0.01 


0.01 


9 


8 


0.1201 


0.0839 


32 


0.01 


0.01 


0.01 


10 


7 


0.1499 


0.0658 


64 


0.05 


0.02 


0.03 


11 


6 


0.1838 


0.0453 


128 


0.17 


0.09 


0.10 


12 


6 


0.2048 


0.0351 


256 


0.72 


0.38 


0.30 


13 


5 


0.2211 


0.0226 


512 


2.93 


1.53 


1.07 


13 


5 


0.2374 


0.0194 


p = 3 


8 


0.00 


0.00 


0.00 


9 


9 


0.1201 


0.1201 


16 


0.01 


0.01 


0.01 


10 


9 


0.1560 


0.1148 


32 


0.02 


0.01 


0.02 


12 


8 


0.1839 


0.0988 


64 


0.10 


0.04 


0.06 


13 


8 


0.2104 


0.0900 


128 


0.41 


0.16 


0.20 


13 


8 


0.2363 


0.0858 


256 


1.76 


0.66 


0.72 


14 


8 


0.2514 


0.0828 


512 


7.45 


2.56 


2.35 


14 


7 


0.2644 


0.0706 


p = 4 


8 


0.03 


0.01 


0.00 


11 


11 


0.1686 


0.1686 


16 


0.01 


0.01 


0.01 


12 


11 


0.2073 


0.1665 


32 


0.05 


0.02 


0.03 


13 


9 


0.2419 


0.1248 


64 


0.24 


0.11 


0.10 


14 


9 


0.2549 


0.1054 


128 


1.07 


0.32 


0.43 


15 


8 


0.2688 


0.0884 


256 


4.47 


1.23 


1.09 


15 


7 


0.2924 


0.0648 


512 


18.79 


5.30 


4.13 


16 


7 


0.3061 


0.0534 



forming the inter-grid transfer operators. It is intended to help the reader in writing optimized/fast 
computer programs. 

(2) The construction of B -spline basis functions at coarse level from the linear combination of fine 
basis functions is provided. For p = 2,3,4, and with C° and C p ~ l continuities, these transfer 
operators (from fine level to coarse level) are given in matrix form for a multilevel mesh. These 
operators can also be used to generate restriction operators in multigrid methods. 

(3) The transfer operators are also provided for NURBS basis functions. The formulation of NURBS 
operators is given in terms of B-spline operators and weights. 

(4) The construction of hierarchical spaces for B-splines (NURBS) is presented. Hierarchical spaces 
are constructed as direct sum of coarse spaces and corresponding hierarchical complementary 
spaces. We have presented matrix form of these operators. As the choice of hierarchical com- 
plementary spaces is not unique, we have provided two different choices of these operators for 
each of C°- and C p_1 -continuity of basis functions. 

(5) For a given polynomial degree p, AMLI cycles are of optimal complexity with respect to the 
mesh refinement. Also, for a given mesh size h, AMLI cycles are (almost) p-independent. We 
provided numerical results for a square domain, quarter annulus (ring), and quarter thick ring. 
The iteration counts, convergence factor, and timings are given for AMLI linear V-, W- and 
nonlinear VF-cycles. Note that, for C p ~ 1 -continuity the linear F-cycle also exhibits optimal 
convergence rates (due to very nice space splitting, which is normally not found in standard 
FEM). The linear and nonlinear AMLI VF-cycle is optimal for all cases except for a particular 
case of degree p = 4 with C° -continuity in three-dimensional problem with first choice of T P ' r . 
For this case, the number of iterations are provided with u = 3 cycle, which is optimal. The 
numerical results are complete for p = 2,3, 4, with C p ~ l and C° continuous basis functions. 

Despite that the condition number of the stiffness matrix grows very rapidly with the polynomial 
degree, these excellent results exhibit the strength and flexibility of AMLI methods. Nevertheless, the 
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Table 6. AMLI methods for Example 2: C° regularity 



with first choice of If' 1 £ 


iven in (59) 


1/h 


tc 


t. 


n it 


P 






N2 


N2 


N2 


p = 2 


8 


0.00 


0.01 


11 


0.1744 


16 


0.01 


0.02 


11 


0.1820 


32 


0.02 


0.05 


11 


0.1791 


64 


0.09 


0.13 


11 


0.1752 


128 


0.40 


0.43 


11 


0.1730 


256 


1.72 


1.52 


11 


0.1717 


512 


7.36 


5.61 


11 


0.1704 


p = 3 


8 


0.00 


0.01 


13 


0.2237 


16 


0.02 


0.04 


14 


0.2507 


32 


0.08 


0.11 


14 


0.2584 


64 


0.34 


0.39 


14 


0.2632 


128 


1.49 


1.43 


14 


0.2649 


256 


6.35 


5.37 


14 


0.2648 


512 


27.51 


20.83 


14 


0.2638 


p = 4 


8 


0.01 


0.03 


22 


0.4319 


16 


0.05 


0.11 


24 


0.4516 


32 


0.22 


0.38 


24 


0.4563 


64 


0.92 


1.34 


24 


0.4591 


128 


4.21 


5.03 


24 


0.4609 


256 


18.28 


19.79 


24 


0.4639 


512 


76.62 


81.78 


25 


0.4644 



with second choice of Tf' 


given in (60) 


1/h 


tc 


ts 


nn 


P 






N2 


N2 


N2 


p = 2 


8 


0.00 


0.01 


10 


0.1445 


16 


0.01 


0.02 


10 


0.1510 


32 


0.04 


0.04 


10 


0.1478 


64 


0.15 


0.15 


10 


0.1463 


128 


0.62 


0.52 


10 


0.1437 


256 


2.65 


1.93 


10 


0.1419 


512 


11.05 


7.72 


10 


0.1401 


p = 3 


8 


0.01 


0.01 


11 


0.1647 


16 


0.02 


0.03 


11 


0.1780 


32 


0.09 


0.09 


11 


0.1845 


64 


0.39 


0.33 


12 


0.1883 


128 


1.63 


1.21 


12 


0.1922 


256 


6.98 


4.52 


12 


0.1938 


512 


28.76 


17.94 


12 


0.1940 


p = 4 


8 


0.01 


0.02 


11 


0.1660 


16 


0.07 


0.05 


11 


0.1758 


32 


0.32 


0.19 


11 


0.1789 


64 


1.39 


0.70 


11 


0.1785 


128 


5.99 


2.64 


11 


0.1774 


256 


25.31 


10.49 


11 


0.1765 


512 


99.22 


43.15 


11 


0.1757 



TABLE 7. AMLI methods for Example 3: First choice of Tf ' r given in (57) with C p 1 regularity 



'fe+i 



1/h 


tc 


ts 


n it 


P 






LI 


N2 


LI 


N2 


LI 


N2 


p = 2 


4 


0.00 


0.00 


0.00 


8 


8 


0.0899 


0.0899 


8 


0.04 


0.01 


0.01 


12 


10 


0.1913 


0.1438 


16 


0.60 


0.10 


0.10 


13 


10 


0.2400 


0.1484 


32 


7.18 


1.09 


0.89 


15 


10 


0.2694 


0.1346 


64 


* 


* 


* 


15 


9 


0.2830 


0.1168 


p = 3 


4 


0.00 


0.00 


0.00 


10 


10 


0.1415 


0.1415 


8 


0.15 


0.02 


0.03 


14 


13 


0.2492 


0.2304 


16 


2.84 


0.27 


0.24 


15 


11 


0.2923 


0.1862 


32 


35.61 


2.79 


2.21 


17 


11 


0.3215 


0.1762 


64 


* 


* 


* 


17 


11 


0.3349 


0.1738 


p = 4 


4 


0.01 


0.01 


0.01 


10 


10 


0.1443 


0.1443 


8 


0.52 


0.06 


0.07 


16 


16 


0.3027 


0.3040 


16 


14.81 


0.82 


0.85 


20 


17 


0.3900 


0.3324 


32 


213.74 


8.82 


7.55 


21 


15 


0.4067 


0.2927 


64 


* 


* 


* 


21 


14 


0.4042 


0.2546 



25 



Table 8. AMLI methods for Example 3: with C° regularity 



with first choice of Tp' given 



in (59) 



1/h 



16 

32 



16 

32 



16 

32 



N2 



n it 



N2(N3) 



V 



0.01 
0.11 
1.25 
12.06 



0.01 
0.05 
0.52 
4.51 



12 (12) 
15(15) 
16(15) 
16 (15) 



P 



0.07 
1.09 

12.23 
114.77 



0.04 
0.50 
5.13 
49.48 



18(18) 

23 (22) 
26 (23) 
28 (23) 



P 



0.39 
5.84 
64.09 



0.45 
4.36 

47.27 



48 (48) 
54 (50) 
64 (51) 
73 (51) 



N2 



0.2124 
0.2904 
0.2996 
0.3022 



0.3527 
0.4408 
0.4919 
0.5164 



0.6770 
0.7081 
0.7497 
0.7764 



with second choice of Tf'' 


given in (60) 


1/h 


tc 


ts 


THt 


P 






N2 


N2 


N2 


p = 2 


4 


0.37 


0.31 


11 


0.1753 


8 


0.32 


0.13 


13 


0.2212 


16 


4.29 


0.76 


13 


0.2250 


32 


33.13 


7.44 


13 


0.2261 


p = 3 


4 


0.09 


0.03 


14 


0.2663 


8 


1.42 


0.34 


16 


0.3092 


16 


15.72 


3.30 


17 


0.3342 


32 


123.05 


32.24 


18 


0.3415 


p = 4 


4 


0.98 


0.23 


16 


0.2987 


8 


13.39 


1.84 


18 


0.3465 


16 


144.03 


17.32 


18 


0.3560 


32 


* 


* 


18 


0.3577 



rigorous local analysis of the CBS constant 7, particularly due to the overlapped support of B-splines, 
is not a straight forward task, and is still an open problem. We intend to address this issue in our future 
work. 
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